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Abstract 

We present a program for the numerical evaluation of scalar integrals and tensor 
form factors entering the calculation of one-loop amplitudes which supports the use 
of complex masses in the loop integrals. The program is built on an earlier version of 
the golem95 library, which performs the reduction to a certain set of basis integrals 
using a formalism where inverse Gram determinants can be avoided. It can be used 
to calculate one-loop amplitudes with arbitrary masses in an algebraic approach as 
well as in the context of unitarity-inspired numerical reconstruction of the integrand. 
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LONG WRITE-UP 



1 Introduction 



Collider experiments at the TeV scale, in particular the LHC experiments, are 
expected to shed light on the mechanism of electroweak symmetry breaking 
and possibly guide us to a more complete theory of fundamental interactions 
than the present Standard Model. In order to achieve these goals, predictions 
for signal as well as background processes should be well under control, neces- 
sitating calculations at next-to-leading order (NLO) accuracy or beyond. 

Over the last few years, enormous progress has been made to push the cal- 
culation of NLO corrections towards a higher number of particles in the final 
states, i.e. to "multi-leg" amplitudes, in QCD as well as in the electroweak 



sector. For reviews see e.g. |5|l6] . 

Nowadays, the efforts are also increasingly being focused on the goals of 
automating multi-leg one-loop calculations and making them publicly avail- 
able. Recent public programs with emphasis on multi-leg one-loop calcula- 
tions are e.g. Cut Tools [7], Samurai Pl9] . other public programs which have 
been optimized for less than four particles in the final state are e.g. Fey- 
nArts/FormCalc pmHITTB] . MCFM fT5], VBFNLO 



An important ingredient for such programs is an integral library containing 
the one-loop integrals which are the basic building blocks of any one-loop am- 
plitude unless it is calculated purely numerically. Several libraries are available 
to date: FF [T1IT5] . Looptools [H], QCDLoop [I6], OneLOop [2], golem95 [IT], 
Hexagon. F |18j. A code for the calculation of one-loop four-point functions 
with complex masses (DOC) can be found in [T^. The latter has been inte- 
grated into the Loop Tools library [5] where the complex version of infrared 
finite integrals with less than four legs are already implemented. A complete 
set of scalar four-point integrals, both in dimensional and in mass regular- 



isation and valid also for complex masses can be found in [2^ in analytic 
form. 

The calculation of scalar one-loop integrals has a long tradition of pioneering 
work, see e.g. |21|22|14|23f24f25|26] . For processes involving unstable particles, 
these integrals are also required for complex internal masses, in order to be able 
to work within the so-called "complex- mass scheme" developed in Refs. [27p8] . 
For calculations beyond one loop, complex values for invariants derived from 
external momenta are also required [29|30] . but we will concentrate on one- loop 
corrections here. 
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In this article, we present an extension of the tensor and scalar library of 
Ref. [17] to integrals with arbitrary masses, in particular also complex masses. 
Furthermore, we extend the approach which was previously based solely on 
form factors as building blocks of the amplitude to an approach which is 
useful in the context of reconstruction of the integrand using D-dimensional 
unitarity. 

This article is organized as follows. In Section 2, we review the theoretical 
background, with particular emphasis on the treatment of potential numeri- 
cal instabilities. In subsection 2.4 we give an example to demonstrate how the 
introduction of complex masses can cure Landau singularities stemming from 
on-shell massive particles in the loop. Subsection 2.5 is dedicated to the new 
feature of golem95 to be used in the context of a numerical reconstruction 
of amplitudes at the integrand level. Section 3 gives a brief overview of the 
software structure, while a detailed description of the individual software com- 
ponents is provided in Section 4. The installation instructions can be found in 
Section 5, followed by a listing of the examples which are new in this version 
in Section 6, before we conclude. 



2 Theoretical background 



The program is an update of the tensor and scalar integral library described in 
more detail in Ref. p2], based on the formalism developed in Refs. [3T]l32] to 
reduce tensor integrals to a convenient set of basis integrals. Similar reduction 
schemes can be found e.g. in Refs. |25|33|34|35|36|37|38j . Here we will describe 
the theoretical framework only briefly and focus on the new features of the 
program. 



2.1 Form Factors 



Tensor integrals can be divided into a part containing the Lorentz structure 
and a part consisting of scalar quantities, which we call form factors. 

We define an A^-point tensor integral of rank r in D = A — 2e dimensions as 



TD,fii...firf^ \_ f d, q ?ai • • • Qar /IN 

yai,...,ar)— / . ,n n , ■s:\ /„2 „„2 , ■s:\ 



i TT^/"^ {ql — ml + i6) . . . {q% — m% + i6) 
where qa = q + ra, qis the loop momentum, and is a combination of external 
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momenta. Using the shift invariant vectors 



ij i j 



we can write 



(2) 



ji,...,jres 

jl,...,jr-2eS 



A: 



'Jr-2- 



jl,---Jr-4e5 



Jr-4,- 



{ai...ar} 



iS). (3) 



The notation [• • stands for the distribution of the r Lorentz indices 

jii, and the momentum labels to the vectors A^'„. and metric tensors in 
all distinguishable ways. S denotes an ordered set of propagator labels, corre- 
sponding to the momenta forming the kinematic matrix S, defined by 



i,j e{l,...,7V}. 



(4) 



The form factors are linear combinations of so-called reduction coefficients 
derived from the matrix S and A'^-point integrals with N < A. The kinematic 
matrix S is related to the Gram matrix Gij — 1, . . . , N — 1 for r^v = 0) 

by 



AT 



detG= (-1)^+^5 detS , B = ^^j^ ■ 



(5) 



2.2 Integrals 



The golem95 program uses the fact that tensor integrals are related to Feyn- 
man parameter integrals with Feynman parameters in the numerator. A scalar 
integral, after Feynman parametrisation, can be written as 



^ ■' i=l 1=1 



-N 



N 



(6) 



5 



The general relation between tensor integrals and parameter integrals with 
Feynman parameters in the numerator is well known [39|25II32] 



j-D,fJ.l...flr 

'at 



N 

E 



[r/2] 



m=0 
{Mi-Mr} 



TD+2m 
'■N 



(j'l • • • ) jr~2m ; S) 



(7) 



where . . . ,jr-2m ; S) is an integral with Feynman parameters in the 

numerator, [r/2] stands for the nearest integer less or equal to r/2 and the 
symbol indicates that m powers of the metric tensor are present. Feynman 
parameter integrals corresponding to diag rams where propagators ji, . . . , jVrt 
are omitted, or pinched, with respect to the "maximal" topology can be defined 
as 



mji, . . . , jr; s\{h,..., im}) = (-i)^r(iv - -) 

r ^ D/2-N 

/ n dz, 5(1 - E ^k) 5izi,) . . . 5{z,Jz,, . . . z,^ (le) . (8) 



The program golem95 reduces the integrals internally to a set of basis in- 
tegrals, i.e. the endpoints of the reduction (they do not form a basis in the 
mathematical sense, as some of them are linearly dependent). The choice of 
the basis integrals can have important effects on the numerical stability in cer- 
tain kinematic regions, as will be explained below. Our reduction endpoints 
are 4-point functions in 6 dimensions J4, which are IR and UV finite, UV 
divergent 4-point functions in D -|- 4 dimensions, and various 2-point and 3- 
point functions, some of the latter with Feynman parameters in the numerator. 
This provides us with a convenient separation of IR and UV divergences, as 
the IR poles are exclusively contained in the triangle functions. Explicitly, our 
reduction basis is given by integrals of the type 
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/f «00 - -r (2 - f ) /; n *(i - g . 
/f-O.) - r (2 - f ) I n -'^^ *(i - E ■ 

(9) 

where r™^^ = 3, as well as I^, I^^'^-, I^^"^-, I^^^ with no Feynman parameters 
in the numerator, and two-point and one-point functions. 

Note that /-f "^"^ and I^^^ are UV divergent, while can be IR divergent. In 
the code, the integrals are represented as arrays containing the coefficients of 
their Laurent expansion in e = (4 — D)/2. 

We would like to emphasize that the program also can be used as a library 
for scalar master integrals. 

2.3 Treatment of potential numerical instabilities 

2.3.1 Spurious singularities due to inverse Gram determinants 

Further reduction of the integrals in eqs. ^ to integrals with no Feynman 
parameters in the numerator introduces factors of 1/-B, i.e. inverse Gram de- 
terminants. A particular feature of golein95 is the fact that the above inte- 
grals are not reduced to scalar basis integrals in cases where B becomes small, 
thus avoiding problems with small inverse determinants. In these cases, the 
above integrals are evaluated numerically. As B = (— 1)^"*"^ det(G')/ det(iS) is 
a dimensionful quantity, the switch to the numerical evaluation of the basis 
integrals is implemented such that the value of the dimensionless parameter 
B is tested, where 

B = B X (largest entry of 5) . (10) 

If 5 > B'^^^, the reduction is performed, else the program switches to the 
direct numerical evaluation of the integral. The default value is B'^^^ = 0.005. 
In particular, we use a certain one-dimensional parameter representation here, 
obtained after performing two integrations analytically. In this way one can 
use deterministic integration routines, leading to a fast and precise numerical 
evaluation. The relative error to be achieved in the numerical integration has 
been set to the default value 10~^. If this precision has not been reached, the 
program will write a message to the file error.txt. 
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This switch to a direct numerical evaluation will be done automatically for all 
triangle integrals, and for box integrals with massless propagators and up to 
three off-shell legs. For box diagrams with massive internal propagators the 
one-dimensional parameter representation is not yet implemented, but will be 
provided in a forthcoming version. 

In some cases, calculating in double precision Fortran may not be sufficient. 
The code is designed such that it can be compiled in quadruple precision as 
well. 

2.3.2 Landau Singularities 

After Feynman parametrisation and momentum integration, the denominator 
V of any one-loop integral is given hy V = —z ■ S ■ z /2 — i5 (see eq. 
Necessary conditions for a Landau singularity to occur in a one-loop integral 
thus can be expressed as 

det5 = 0, 2i>0, z-S-z <Q . (11) 

The leading Landau singularity, corresponding to det 5 = and -Zj > for all z, 
occurs if all particles in the loop go simultaneously on-shell. Sub-leading Lan- 
dau singularities occur if a sub-matrix of S has a vanishing determinant and 
at least one of the Feynman parameters Zi is zero, corresponding to pinched 
diagrams. Leading Landau singularities are also called anomalous thresholds 
in the literature |10|^ . this term stemming from the fact that they only oc- 
cur in the physically allowed phase space region if a number of rather special 
kinematical conditions are fulfilled. An example will be given below. 

2.4 Complex masses 

In order to demonstrate how the introduction of complex masses, standing 
for a width in the case of unstable particles in the loop, regulates the Landau 
singularities, we choose the example of a box diagram contributing to the 
production of a heavy neutral Higgs boson and a hb pair in gluon fusion in the 
context of supersymmetry, where the loop contains two squarks (sbottoms) 
and two neutralinos, as shown in Fig.[Tl 

Denoting the momentum of the Higgs boson by and the momenta of the 
two b-quarks by ^3 and p4, the kinematic matrix iS*^^-* associated with the 
four-point diagram is symmetric under S35 -H- S45 {sij = {pi +pj)'^). To make 
the singularity structure apparent, we will scan the different thresholds and 
singularities as a function of the invariant S45. The 6-quarks are taken to be 
massless in this example. Similar investigations have been worked out in detail 
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Fig. 1 . A diagram where all propagators can go simultaneously on-shell to develop a 
leading Landau singularity which can be regulated by introducing complex masses. 

in piiH] . 

Solving the phase space constraint det G > {Gij = 2piPj,i = 3, 4, 5) for S45 
leads to the following boundaries for the physically allowed phase space: 



Ml — < S45 <Ml + si2 - S35 
S35 



:i2) 
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The leading Landau singularity is characterised by all particles in the loop 
going simultaneously on-shell, leading to det S^'^^ = 0. The determinant of the 
kinematic matrix S^j'^ can be written as [l3j 



det = A(S35, ml Ml) (545 - 4)' + det Sf^ det S'g (13) 
A(x, z) = + + — 2 {xy + xz + yz) . 

The determinants det Sf^ and det Sf^ correspond to diagrams where the prop- 
agators r2 and respectively in Fig.[T] are pinched. The Kaellen function 
A(s35, m?,m^) is (minus) the determinant of the kinematic matrix associated 
with a two-point function where both ri and are pinched. S45 is the solution 
of the equation det 5^'^^ — det Sf^ det Sf^ = 0. Choosing the numerical values 
tuh = 450 GeV, = 800 GeV, = 200 GeV, = 1700 GeV, fixing S35 to 
S35 = 2(m? -|- m^) and combining with the phase space constraints of eq. (fT2l) 
and the requirement z ■ S ■ z < 0, we encounter the following discontinuities, 
shown in Fig.[2l 

• a normal threshold at ^545 = (m^ + m^) = 1000 GeV, where A(s45, m?, M^) 
vanishes, corresponding to the production of a squark and a neutralino in 
the cut two-point diagram associated with a pinch of both r2 and r4 in the 
box diagram. 



9 



• an anomalous threshold at ^/s^ — 1012.7 GeV, corresponding to 
det^g) = 0. 

• an anomalous threshold at ^/s^ — 1038.1 GeV, corresponding to 
detSf} = 0. 

'4 

• a leading Landau singularity at ^/s45 ~ 1078.4 GeV, corresponding to 
det5(^) = 0. 




Fig. 2. Singularity structure of the scalar four-point function A^''^ (real masses) con- 
tained in the diagram of Fig.[T] for itlh = 450 GeV, rriq = 800 GeV, = 200 GeV, 
= 1700 GeV, 900 GeV < S45 < 1200 GeV. 

Introducing complex masses moves the poles due to propagators going on-shell 
away from the real axis and therefore regulates the leading Landau singularity, 
as can be seen from Fig. [31 

We should remark that Landau singularities due to massive (unstable) par- 
ticles in the loop going on-shell usually behave like l/V det S as det5 — )■ 0. 
This is in contrast to cases with several massless particles, where det S usually 
has several zero eigenvalues. An interesting case is the six-photon amplitude, 
where not only det 5, but also its derivative with respect to the invariants 
involved vanishes at the singular point, meaning that the singularity is not 
integrable anymore. Introducing an imaginary part as advocated above would 
certainly not help in this case, but in the six-photon example we are saved 
by the gauge structure, which leads to numerators taming the singularity 
structure when the individual contributions are combined to physical helicity 
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amplitudes. This has been worked out in detail in Refs. |HP5f5] . 
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Fig. 3. Singularity structure of the scalar four-point function (real masses) contained 
in the diagram of Fig.[T]for mn = 450 GeV, m? — )• m? — i rriqTq, — )■ — i m-x^^, 



Tq = 3.5 GeV, = 1.5 GeV 



2. 5 Tensorial Reconstruction of the Integrand 



The library golem95 in its original version[T7], when used in amplitude cal- 
culations, relies on some user generated code expressing a diagrani ^ I in terms 
of form factors (A^'*", i?^'^', C^'*') as defined in Eq. ([3]). This requirement of 
the previous version of golein95 restricted the applicability of the library to 
algebraic methods for generating the amplitude. In order to be able to use 
the library in the context of a numerical reconstruction of the integrand, we 
included new features described in the following. 

The general structure of a one-loop diagram can be written as 

d^g J\f{q) 



(^q2 _ ^2 _|_ . . . (^q2^ _ _|_ ^^-j 



(14) 



^ Actually, it is not required to organize the calculation in terms of diagrams. The 
term diagram here refers to a set of terms of an amplitude sharing some common 
loop propagators. 
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where qa = q + Ta- After decomposing the D- dimensional vector q into its 
projection q onto the physical 4-dimensional Minkowski space and the radial 
component /i of the {D — 4)-dimensional orthogonal space, such that q^ = 
— /i^, we can express the numerator function Af{q) in terms of the following 
tensor structure, 



R-2a 



(15) 



0=0 



r=l 



The coefficients C'^'^ can be determined numerically by the algorithm described 
in Ref. pj. The contraction of the coefficients with the tensor integrals can be 
carried out in a process independent and numerical way. Therefore this tenso- 
rial reconstruction of the integrand provides a way of processing both diagrams 
from algebraic constructions and purely numerical input, using the numerator 
function A^(g, /i^) and the set of denominators as the only common source of 
information. We have implemented the reconstruction of the coefficients and 
their contraction with the tensor integrals as a part of golem95. The new 
interface is described in Section 14.21 

The contraction with the tensor integrals also requires the implementation of 
integrals with explicit /x^-dependence, which can be computed from the poles 
of known form factors. 



d^g 



i7r-D/2 ^q_^ _ ^2 _|_ ^(J-) . . . _ ^2^ _|_ 

(A-D) 




{ai...a,_2} Ji--''-^^ ^ 



+ E 

jl,...,jr-4eS 



(16) 

{ai...ar-2} Ji- Jr t I 



In practice, the above equation simplifies greatly since the only non-zero con- 
tributions come from the form factors B^'"^, B^'"^, B^'^ and C^'^. The only non- 
zero integral with /x^ in the numerator is the box 



d^g 



i7!-D/2 _ ^2 _|_ ^(j-) . . . 



ml + iS) 

{D-A){D-2)C^'\S) + 0{D-A). (17) 



3 Overview of the software structure 



The structure of the golem95 program is the following: There are four main 
directories: 
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(1) src: the source files of tlie program 

(2) demos: some programs for demonstration 

(3) doc: documentation which has been created with robodoc [IB] 

(4) test: supplements the demonstration programs, containing files to pro- 
duce form factors with user-defined kinematics. The user can specify the 
rank, numerator, numerical point etc. via a steering file. 

The subdirectory structure is the same as described in [T7] . 



4 Description of the individual software components 

We focus here on the new features, for more details on the software components 
which are the same as in version 1.0, we refer to [17j and to the documentation 
contained in the program. 

4-1 Form factor evaluation 

A typical setup for the calculation of form factors is: 

call initgolem95(dim) 
... fill matrix S . . . 
call preparesmatrixO 
. . . evaluate form factors . . . 
call exitgolem95() 

The three subroutines act as interfaces to the library. They are used to initial- 
ize and manipulate objects and features of the program. They are implemented 
in the module matrice_s and will be described in more detail in the following. 
To pass the desired masses and momenta, the user assigns values to the object 
s_mat, which are the entries of the kinematic matrix S given in eq. (jl]). The 
form factors are implemented in the respective modules in src/form_f actor. 
We will sketch the internal structure of the library and the distinction between 
real and complex masses in the following. 

4-1.1 The module matrice_s 

This module is located in src/kinematic/matrice_s . f 90. It is used to reserve 
and free memory for the kinematic matrix S and related objects, as well as for 
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the computation of the inverse matrices. The three macro functions described 
above are: 

subroutine initgolem95(dim, opt_set) allocates memory for two arrays 
s_mat_c and sjnat_r of dimension dim and rank two. These arrays repre- 
sent internally the complex valued matrix S and its real part. Additional 
memory is reserved for the inverse of S and its submatrices as well as for 
parameters appearing in the reductions of the form factors. A public pointer 
s_mat is associated with s_mat_c. This object is used to assign the matrix 
iS, both in the real mass and complex mass case. An additional argument 
opt_set can be given, which is an integer array reflecting the numbering of 
the propagators. The default value is an integer vector of range 1 to dim. 
This subroutine also initializes the caching system. 

subroutine preparesmatrixO creates the derived type object s_mat_p. This 
consists of two pointers to s_mat_c and s_mat_r as well as two integers encod- 
ing information about the entries of complex masses and vanishing masses. 
We describe the usage of this object in more detail in section 14.1.21 The 
arrays used later for the inverse matrices and sub-matrices are calculated 
in this routine. If the user defines values for the matrix S which are purely 
real, only the real arrays are needed and calculated. 

subroutine exitgolem95() deallocates all arrays, nullifies all pointers and 
clears the cache. 

Note that the subroutines initgolem95 and exitgolem95 also have been de- 
signed to make the setup more user-friendly, so the calls to the subroutines 
allocation_s, allocate_cache, init_invs, deallocation_s, clear_cache 
and the definition of set_ref in version 1.0.0 are now obsolete. 

4.1.2 Internal structure 

Internally, the library is subdivided into three layers, implementing the reduc- 
tion from general form factors down to specific integrals. The global derived 
type object s_mat_p is passed from one layer to the next. The information 
encoded in the integer bits of s_mat_p, is used to direct the evaluation of the 
form factors to the specific integrals needed. 

Upon a form factor call, reductions are performed in the respective modules 
located in src/form_f actor. These reductions only need information from the 
inverse matrices and related objects created with preparesmatrixO. 

From the form factor modules, generic iV-point functions are called, passing 
the derived type object s_mat_p. At this level, further reductions might be 
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performed. Up to this point, the only information needed from the kinematic 
matrix S is the position of the zero-mass entries, which is encoded in the 
integer sjnat_p°/ob_zero. In the current implementation, we choose massless 
six- dimensional boxes and massive four-dimensional boxes as our basis of the 
reduction. A comparison of b_zero with the pinched propagators gives a fast 
case distinction. 

Only in the last step, when the generic A^-point functions call specific integrals, 
does the (sub-) matrix S need to be passed over. If the latter contains complex 
masses, which again can be quickly determined by comparing s_mat_py„b_cmplx 
with the set of remaining propagators, an implementation of the integral for 
complex masses is called with s jiat_p°/„pt_cmplx, otherwise the version for real 
masses is called with s_mat_p%pt_real. This structure of the library allows for 
an efficient reduction of the form factors to the basis integrals, without a full 
knowledge of the kinematic matrix S. 

4-2 Interface for tensorial reconstruction 

The tensorial reconstruction described in Section 12.51 has been implemented 
in two modules, tens_rec and tens_conib. In the simplest case the user only 
needs to access the function evaluate from the module tens_comb. For more 
advanced applications of the interface the user will have to use functions and 
data types of both modules. 

4.2.1 The module tens_rec 

The module tens_rec is located in src/interf ace/tens_rec . f 90 and con- 
tains all routines required for the reconstruction of the tensor coefficients and 
for the evaluation of reconstructed numerators. 



4.2.1.1 Data types The module tens_rec defines the data types 

coef f _type_l, . . . , coef f _type_6, 

with components ck {0 < k < min(4, R)) to store the coefficients of a numera- 
tor of maximum rank R. In particular, the component ck is a two dimensional 
complex array which stores the coefficients of the monomials with k distinct 
components of q. The first index of the array labels the (J^ ways of choosing 

k of the four components of q. The second index labels the J2^=i (lli) — [T) 
monomials of the corresponding polynomial. A certain order of the array en- 
tries should not be assumed by the user. 
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4.2.1.2 Functions and subroutines In all functions and subroutines 
which have a parameter numeval it should point to a function with the fol- 
lowing signature: 

interface 

function numeval (Q, mu2) 

use precision_golem, only: ki 
real(ki), dimension(0 : 3) , intent(in) :: Q 
real(ki), intent(in) :: mu2 
complex (ki) : : numeval 
end function 
end interface 

For all momenta the order k(0 : 3) = {E, x, y, z) is assumed. 

subroutine reconstructr (numeval, c0,cl,c2) forr = 1, 6 reconstructs 
the tensor coefficients of the numerator function given by numeval. The pa- 
rameters cO, cl and c2 are output parameters; cl and c2 are optional 
and only available in the subroutines where r > 2. The argument cO is of 
type coeff _type_r and contains the constant part (with respect to /i^, i.e. 
a = in eq. f|T5l) ) of the polynomial. The arguments cl and c2 are of type 
coeff _type_(r — 2) and contain the coefficients of the /i^ and /i"^ part of the 
polynomiaO]- 

subroutine print_coeffs (coeff s , unit) is overloaded and can take any 
derived coefficient type as its first argument. The second argument is op- 
tional and has the default value unit=6 (standard output). This routine 
pretty-prints the coefficients to the given file or device. 

pure function tensevalr(Q, coeff s, max_k) evaluates the polynomial given 
by the coefficients coeff s for the given real momentum Q. The optional ar- 
gument max_k is only used internally and should not be assigned by the user. 

pure function ctensevalr(Q, coeffs) evaluates the polynomial given by 
the coefficients coeffs for the given complex momentum Q. 

The subroutines reconstructr and ctensevalr can together serve as an im- 
plementation of a presampling for programs of unitarity based reductions at 
the integrand level [7|8 ] or as an implementation of a rescue system when 
the unitarity based method fails due to vanishing Gram determinants. Both 
options are described in more detail in [H]. 

4-2.2 The module tens_comb 

The module tens_comb is located in src/interf ace/tens_comb . f 90 and con- 
tains all routines for the contraction of reconstructed tensor coefficients with 

^ There is no coeff _type_0: complex (ki) is used instead. 
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tensor integrals. It also contains the convenience function evaluate which 
combines the steps required for reconstruction and contraction. 



4.2.2.1 Functions and subroutines 

function evaluate (numeval, momenta, set, rank) combines tensorial re- 
construction and the contraction of the coefficients with the tensor integrals. 
The argument momenta is of dimension( : ,0:3) and contains the momenta 
ri, . . . , ttv. The argument set can either be an integer array or an integer 
number (generated by the function packb) and denotes the set of pinched 
propagators; a value of zero can be used if no propagators are pinched. It 
should be noted that the argument momenta also includes the momenta be- 
longing to pinched propagators. The last argument, rank, is optional. It 
specifies the rank of the numerator; if omitted it is assumed that rank = A^. 

function contractiV_r(coeff s, momenta, set) contracts the coefficients 
coeff s of type coeff _type_r with the A^-point rank-r tensor integral. The 
arguments momenta and set are defined as in evaluate with the restriction 
that set has to be an integer number. 

function contractA^_rsl (coeff s , momenta, set) contracts the coefficients 
coeff s of type coeff _type_(r — 2) with the corresponding A^-point tensor 
integral containing fi^ in the numerator. The arguments momenta and set 
are defined as in evaluate with the restriction that set has to be an integer 
number. 

function contractA^_rs2 (coeff s , momenta, set) contracts the coefficients 
coeff s of type coeff _type_(r — 2) with the corresponding A^-point tensor 
integral containing //^ in the numerator. The arguments momenta and set 
are defined as in evaluate with the restriction that set has to be an integer 
number. 

The function evaluate determines A^ from the size of momenta and the number 
of pinches in set. It then combines the according calls to reconstructr and 
contract A^_r, contract A^_r si and contract A^_rs2. The functions evaluate 
and contract . . . require that the matrix s_mat and the cache are set up 
properly by the sequence which one would have to call when evaluating form 
factors. 



5 Installation instructions 

The program can be downloaded as an archive golem95-1.2.0.tar.gz from the 
following URL: http://projects.hepforge.org/golem/95/. The installa- 
tion instructions given below also can be found in the Readme file which comes 
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with the code. Information and updates of the program can also be found at 
http : / /pro j ects . hepf orge . org/golem/trac/wiki/golem 95C[ 

The installation setup is based on autotools [47j. To install the golem95 li- 
brary, type the following commands: 

./configure [ — pref ix=mypath] [ — precision=quadruple] [FC=compiler] 

[F77=f ortran77compiler] 

make 

make install 

The — prefix option denotes the installation prefix, under which the direc- 
tories lib/ and include/ are generated. If no option is given, on a Linux 
system the configure script would choose — pref ix=/usr/local. The argu- 
ment — precision selects if double or quadruple precision should be used in 
the library; it should be noted that quadruple precision is not supported by all 
Fortran compilers and therefore — precision=double is the default value. If 
the variable FC is not set the first fortran compiler which is automatically de- 
tected will be used. Another variable commonly used is FCFLAGS which allows 
one to pass compiler flags to the Fortran compiler. The fortran77 compiler is 
only needed to run the demo file demos/demo_LT. f . 

In addition, it is possible to call the finite scalar box and triangle integrals 
with internal masses from LoopTools. The setup to do so is automated, all the 
user has to do for this option is (a) install LoopTools, and (b) use the option 
[ — with-looptools=your_path_to_libooptools . a] for the configure script, 
i.e. type the following commands: 

. / configure [ — pref ix=mypath] [ — with-looptools=your_path_to_libooptools . a] 

[ — precision=quadruple] [FC=compiler] [F77= f ortran77compiler] 

make 

make install. 



6 Examples 

Examples can be found in the subdirectory demos. The program demos/ 
demo_cmplx_masses . f 90 shows how the library is used to evaluate form fac- 
tors. For a given matrix S with real or complex masses, six-point form factors 
from rank to rank 6, as well as one- to five-point form factors are evaluated. 
The results are written to the files test_ff (6)_c/r .txt. 

The program demos/demo_tens_rec .f 90 together with the module demos/ 
demo_tens_mod.f 90 demonstrates the use of the function evaluate. In this 
example, a toy amplitude is given by a six-point integral where the numera- 
tor consists of three propagators which are also present in the denominator. 
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The scalar three-point function resulting from the direct cancellation of three 
propagators is compared to the results obtained by expanding the numerator 
into contracted loop and external momenta, leading to rank 6 hexagons if no 
propagators are cancelled, rank 4 pentagons if one propagator is cancelled, 
rank 2 boxes if two propagators are cancelled. 

The program used to scan the Landau singularity as in Fig. [3] can be found in 
demos/SusyLandau . f 90. 



7 Conclusions 

We have presented a program for the numerical evaluation of scalar integrals 
and tensor form factors entering the calculation of one-loop amplitudes, which 
is able to provide results for real as well as complex masses in the loop integrals. 
The program is an extension of on an earlier version of the golem95 library, but 
now also can be used in the context of a unitarity-inspired numerical recon- 
struction of the integrand at the tensorial level. Improvements in the caching 
system and in the user interface also have been made. The program, available 
at http://projects.hepforge.org/goleni/95/, provides a complete library 
of scalar and tensor integrals up to rank six 6-point functions, including the 
option of complex masses. 
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